The present invention relates to a method of analysing images of a deformable object undergoing non-rigid motion. In particular it relates to a method of analysing the image so that desired image features can be detected and tracked through the sequence, and so that the motion of the features can be automatically quantified and analysed.
Basic techniques for analysing images of objects in motion are relatively straightforward for rigid motion, i.e. where the object does not itself deform. However, the analysis of non-rigid motion, where the object deforms in time, is more difficult. Such non-rigid motion occurs in many situations, but a typical one is in the medical or vetinary imaging field where organs of the human or animal body are imaged in real-time. As well as the problem created by the non-rigid motion of the organ being analysed, these imaging applications in particular have the additional problem that the images are very noisy. It is often difficult even for trained operators to find desired image features in the image, and thus reliable automatic detection of the features presents considerable difficulty.
In the field of cardiac imaging, various techniques are used such as multi-gated acquisition scanning (MUGA), fast computed tomography (CT), positron emission tomograph (PET), magnetic resonance imaging (MRI) and echocardiography (i.e. ultrasound imaging). Of these echocardiography is the most widely used because the imaging equipment is relatively cheap to obtain and maintain and the equipment is relatively portable.
In assessing cardiac function the performance of the left or right ventricle of the heart is particularly significant and there has been an increasing interest in obtaining ventricular measurements, such as the chamber dimensions, area, volume and ejection fraction. To provide a more accurate picture of the ventricular function, and in particular to enable assessment of abnormalities which occur in only parts of the ventricular wall, two particular aspects of the motion of the ventricle have proved significant. These are endocardial wall motion, also referred to as wall excursion, and myocardial wall thickening. It has been determined that when the heart muscle becomes ischemic (i.e. deficient in blood), its motion is altered almost immediately. Because abnormalities can be confined to particular parts of the ventricular wall, a systematic method for the assessment of wall motion involves the segmentation of the surface into a number of different segments. The wall can be segmented in various ways, but a useful method is the sixteen-segment anatomical model of the heart proposed by the American Society of Echocardiography and illustrated in FIG. 8 of the accompanying drawings. This is useful in assessing the images derived from two-dimensional echocardiography.
FIG. 8(d) shows an example of a view which is used for analysis in the techniques described below. In assessing cardiac function for instance of the left ventricle clinicians examine the motion and thickening of each segment and try to assess visually the motion of each segment through the heart cycle. One scoring scheme requires the clinician to score each segment as follows:
ScoreGradingCharacterized by1NormalA uniform increase in wall excursion and thickening2HypokineticA reduced (<5 mm) inward systolic wall motion3AkineticAn absence of inward motion and thickening4DyskineticSystolic thinning and outward systolic wall motion
However, this scoring scheme is highly subjective. Thus clinical reporting of echocardiography examination is highly operator-dependent and basically qualitative. Also each segment must be classified as normal or abnormal, so it is difficult for a clinician to indicate within the scoring system subtleties such as only part of a segment being abnormal.
While there is therefore a clear need for an automatic method of detecting and quantifying wall motion and wall thickening, the images are extremely difficult to assess automatically. In the accompanying drawings FIGS. 1(a) to (d) show a typical set of echocardiographic images. FIG. 1(a) shows an image digitized from a video recording; FIG. 1(b) shows an image obtained from one of the latest ultrasound machines; FIG. 1(c) shows a stress echocardiography image of a patient at rest and FIG. 1(d) shows the same patient as in FIG. 1(c) at a peak dose of dobutamine (a drug which mimics the effects of exercise). It will be appreciated that identifying the desired regions of the ventricle is difficult for a human, and that automatic analysis is even more difficult.
Automatic boundary detection of regions in an ultrasound image is available on certain machines manufactured by Hewlett-Packard by a technique known as acoustic quantification (AQ). This technique discriminates boundaries prior to image formation by using discontinuities in the signal returning from the tissue. Pixels with an intensity gradient above a user-defined threshold are marked as boundary points. Pixels labelled as boundary points are then joined together to form connected boundaries. However, FIGS. 2(a)–(d) show that this technique is not always useful. FIGS. 2(a) and (c) show the basic echocardiographic image, and FIGS. 2(b) and 2(d) show the corresponding respective AQ images at different thresholds. It can be seen that the illustrated boundaries do not help assessment of the image at all because they do not accurately follow the real boundaries.
A better method for detecting the inner boundary of the left ventricle (the endocardium) is proposed in the paper “Evaluating A Robust Contour Tracker On Echocardiographic Sequences” by Jacob, Noble, Mulet-Parada and Blake, published in Medical Image Analysis (1997/8 volume 3, number 1, pp 63–75) which is hereby incorporated by reference. As proposed there the inner boundary, endocardium, is modelled by a non-rigid contour (a B-spline) and the variation in the shape of this contour through the echocardiographic sequence (i.e. as the heart contracts and expands) is represented by using a shape-space. This means that the position of the endocardial wall in each image is regarded as being composed of a time varying departure from a defined position e.g. the initial position, the departure being characterised as a time-varying weighted sum of certain basic types of motion of the contour. For instance, a very simple shape-space would characterise the motion of an object as consisting of a certain proportion of rotation and a certain proportion of translation compared to a defined position. Then the only thing which varies with time is the relative amount of the rotation and translation. In analysing echocardiograms a more complicated shape space has been found to be necessary. The paper referred to above uses a principal component analysis (PCA) of the motion of the endocardial wall to find a set of define motions which can efficiently be used as components to approximate the actual motion of the endocardial wall. Again, the only thing which varies through the sequence is the relative weight of the different define motions in each image.
In this technique for detecting and tracking the endocardial wall, the clinician is required first to examine the frames of the image sequence and manually to locate and trace in a few of the frames the endocardial wall. For instance, in a sequence of 60–80 frames the clinician could manually “draw around” the endocardial boundary every fifth frame. A B-spline curve is then fitted to the manually traced contours to provide an approximation of them and a principal component analysis is performed to find the define components of the motion of the contour through the image sequence. Then the whole sequence is reprocessed so that starting from a predefined initial position the position of the endocardial wall in each frame is predicted based on the position in the preceding two frames and the PCA results. The prediction is corrected in each frame by searching for image features (such as intensity edges) representing the actual position of the endocardial wall. When this process is complete, the B-spline curve for each frame can be displayed overlying the image on that frame so that when the sequence is displayed the contour appears to track the endocardial wall through the sequence.
Illustrating this in more detail, it will be recalled that the shape of a B-spline curve is determined by the position of its control points. Thus the movement of a spline curve fitted to the endocardial wall through the image sequence can be expressed entirely as a change from frame to frame of the position of the control points of the spline curve. The x and y coordinates of the control points are conventionally written in a matrix known as a spline-vector Q and as discussed above, the position of the control points in any frame of the sequence can be expressed as an offset from a defined position Q0. The offset, which is time-varying, can conveniently be separated into a time-varying part known as the shape-space vector X and a part representing the type of allowed motions (the main components of the motion), known as the shape matrix W (normally assumed to be constant).
Thus, in matrix notation:—Q=Q0+WX 
In order to find the spline curve which fits to the endocardial boundary in every frame of the sequence (which amounts to finding the position of the control points of the curve in every frame of the sequence) the first step is that the clinician manually draws around the boundary in several frames, for instance every fifth frame. Then a B-spline curve is fitted to the drawn boundary using a user-defined number of control points. FIG. 3(a) shows a schematic representation of a quadratic approximating B-spline with 24 control points used to model the endocardial boundary. FIG. 3(b) illustrates the curve superimposed on a frame of an ultrasound image. A principal component analysis is then performed on the positions of the control points in each of the frames segmented by the clinician to calculate Q0 and W. The aim then is to find the position of the endocardial boundary in all of the frames of the sequence automatically, i.e. without requiring the clinician manually to draw around them. FIG. 4 illustrates the process schematically. The process involves predicting the position of the boundary each frame based on the boundary in the preceding frames. In other words the value of X (the shape-vector) which is the only time varying part is predicted based on the value in the preceding two frames. Then a search is performed around the predicted position to find image features representative of the actual position of the endocardial boundary. In this technique the searches are performed along a plurality of normals spaced along the predicted curve and the image features are identified through known image processing operations, such as looking at the intensity variation along the search line. When the image features corresponding to the boundary have been found the predicted position can be updated and the actual position of the contour (expressed through the position of the control points of the B-spline) is established.
It will be understood that having represented the endocardial boundary as a B-spline, the only time varying part through the sequence is the position of the control points and, in the shape-space representation, the shape-vector X. It will be recalled that the elements of X are the weights of the different types of motion (i.e. the different principal components of the motion) found in the principal component analysis.
FIG. 5 illustrates an example of a principal component analysis performed on four cardiac cycles of an ultrasound image sequence using a B-spline with 14 control points. The six most dominant modes are shown, each is shown as an initial template (thick solid line) and the change in shape represented by that component of the motion is indicated by the thin solid lines. The diagram in the top left of FIG. 5 is the dominant mode and that in the bottom mode is the least dominant. The deformation represented by each mode is shown in an alternative way in FIG. 6 where a flow vector is centred at the start of each span of the spline curve and shows the movement of that part of the curve. Thus the motion of the contour (spline curve) through the image sequence which was analysed can be expressed as a sum of these motions. The position of the curve in any frame represents a certain weighted sum of these motions. The weights are the values of the components of X and thus the motion can be expressed entirely by looking at the time variation of the components of X. The variation of X with time is illustrated for an echocardiogram sequence in FIG. 7. FIG. 7(a) shows the PCA based components for this sequence and FIG. 7(b) shows the values of the weights versus time of each of those components. It is, however, difficult to interpret clinically the significance of these components and weights. For example, the first deformation mode in FIG. 7(a) (top left) appears to be a scaling of the inferior part of the left ventricular boundary. The corresponding plot of the weight of that component illustrates that this motion is basically periodic. But because this component of the motion affects more than just a single part of the boundary (all parts of the boundary move) it does not give a good idea of how any particular region of the wall is moving. Also, some of the information about the motion of that part of the boundary is “encoded” in the other components.
Thus although the principal component analysis, which gives the component magnitudes of the shape-space vector X is very useful in tracking, it does not provide a good basis for automatic interpretation of the results.
It was mentioned above that wall-thickening, known as myocardial thickening is also a clinically significant factor in assessing the condition of the heart. As the heart is beating the ventricle expands and contracts, predominately by periodic thickening of the ventricular wall. If the wall fails to thicken then the ventricular volume will not change by so great an amount and the pumping of blood will be reduced. Thus it would be advantageous to be able to quantitatively analyse the degree of thickening of the ventricular wall. It may be thought that this could straightforwardly be done by detecting the outer (epicardial) boundary of the ventricular wall, in just the same way as the inner (endocardial) boundary is detected above. However, the endocardial boundary (the inner boundary) is a boundary between muscle and blood which have quite different acoustic impedances. In general this means that the endocardial boundary shows up well on an echocardiogram. The epicardial boundary on the other hand is a tissue—tissue interface and so it is very difficult to trace on the image.
Thus even having tracked the endocardial boundary, it is difficult to detect and quantitatively analyse the movement of the epicardial boundary.
The present invention provides techniques which are useful in solving these two problems. Although illustrated in use in analysing echocardiograms of the left ventricle the techniques are not limited to this. They are applicable to analysis of non-rigid motion in two or three dimensions of other deformable objections. Thus they have other medical and vetinary applications as well as being applicable to imaging deformable objects in general and are also applicable to image modalities other than ultrasound.
The first aspect of the present invention provides a method of analysing a sequence of images of an internal body organ in non-rigid motion, comprising the steps of: detecting the boundary of the organ in each image of the sequence; and automatically calculating the amount of movement through the sequence of each of a plurality of clinically significant segments of the detected boundary.
The amount of movement of each of the clinically significant segments, which can be the segments illustrated in FIG. 8, preferably those in FIG. 8(d), can be displayed graphically, for instance as a graph. Further, an average of the amount of movement of that segment can be calculated as a single number representative of the amount of movement of that segment. It is also possible to calculate the variation in the amount of movement in a segment, the greater the variation, the more likely it is that only a part of that segment is normal. It is also possible to calculate and output the maximal excursion of the detected boundary during the motion, for each segment.
Preferably the boundary is detected and tracked by the technique of principal component analysis and fitting of a spline curve as described above.
The amount of movement of the segments can conveniently be found by calculating and outputing for each segment a measure of the amount of movement of the control points controlling the curve within that segment. This measure may be a simple average, or can be weighted in favour of the control points in the middle of each segment.
The variation in the amount of movement within the segment is conveniently found by comparing the amount of movement of the different spline curve control points for that segment.
These measures can easily be obtained from the position of the control points in each frame of the sequence by defining a new shape-space, different from that used in the tracking process, and calculating from the control points the shape-vector corresponding to the different shape-space. The new shape-space can be selected to ensure that each component of the shape-vector represents the amount of movement of control points in a single clinically significant segment only. Then displaying graphically the time varying components of the new shape-vector gives a good indication of the motion of that segment. This aspect of the invention is particularly applicable to analysing ultrasound images of a heart, e.g. of the left or right ventricle.
The invention also contemplates the interpretation of a moving spline curve tracked in one shape-space by using a different shape-space. Thus another aspect of the invention provides A method of analysing a sequence of images of a deformable object in non-rigid motion comprising the steps of detecting a boundary of the object in each of a plurality of frames of the sequence, fitting a spline curve to the boundary in constructing a shape space representation of the movement of spline curve using a first shape space so that the spline curve tracks the boundary, and decomposing the tracking spline curve using a second different, shape space. The different shape-space can be chosen to select a particular attribute of the motion. In other words, a motion tracked using one shape-space need not be interpreted in the same shape-space: a different one—an interpretational shape-space can be used.
Another aspect of the invention involves the modelling of the configuration of the wall of an object having two boundaries as seen in a sequence of images by developing a model of the distance between the two boundaries. Thus rather than modelling each boundary separately, a model is constructed for the distance between them. Thus this aspect of the invention provides a method of analysing a sequence of images of a deformable object in non-rigid motion comprising detecting first and second boundaries of the object in a plurality of frames of the sequence and constructing a shape space representation of the variation through the sequence of the distance between the two boundaries.
The model can be a shape-space of the change in distance between the boundaries, and can be based on a principal component analysis of the way that distance changes. The model can be improved by searching the images to find image features representative of the outer boundary. The model avoids incorrect results such as the outer boundary crossing inside the inner boundary.
Another aspect of the invention provides a method of analysing a sequence of images of a deformable object in non-rigid motion to detect inner and outer boundaries of a wall of the object, comprising the steps of: detecting the inner boundary; and searching outside the inner boundary for image features representing the outer boundary.
Thus because it is known that the inner boundary will be inside the outer boundary, this provides a useful start point for the search of the image features representing the outer boundary.
Preferably a spline curve is fitted to the detected image features representing the outer boundary, e.g. by: manually locating the inner and outer boundaries in only some images of the sequence, calculating a shape-space for the change through the sequence of the distance between the two boundaries, detecting the inner boundary and performing said search outside the inner boundary for image features representing the outer boundary in other images of the sequence; and fitting a spline curve to the detected image features in said other images of the sequence by using said shape-space.
Thus in this method a shape-space representing the distance between the two boundaries is obtained. The use of this shape-space helps to ensure that the calculated outer boundary always lies outside the inner boundary. The distance between the two boundaries in the manually treated frames can be subjected to a principal component analysis which is used as a basis for the shape space. This then provides a convenient model of the deformation of the object wall.
The search for the image features representing the outer boundary can be, for instance, by analysing changes in the image intensity, such as a maximum in the intensity, along search lines projected outwards from the inner boundary. In certain images the plot of the intensity can be extremely noisy and it can be rather difficult to detect a clear maximum. In this case a wavelet decomposition of the profile of image intensity can be performed to smooth the profile.
To obtain a better fit to the actual outer boundary, it is possible to apply extra conditions to the fitting. For instance, detected image features can be down weighted if they imply a curvature of the outer boundary which is too high. Similarly, features can be weighted down if they imply a difference between the inner and outer boundaries which lies outside the shape-space.
This technique is particularly useful for the analysis of ultrasound images, in particular of the heart, e.g. of the left or right ventricle. In that case the distance between the inner and outer boundaries represents the myocardial thickness and the change in that thickness through an image sequence is indicative of the condition of the heart.
Again, the ventricular wall can be segmented according to the FIG. 8 model and the degree of thickening for each separate segment can calculated and graphically displayed, as can the variation within each segment.
It will be appreciated that the methods of the invention are conveniently embodied in a computer program, and thus the invention provides a computer program and computer system for performing the methods above.